knitr::opts_chunk$set(echo = TRUE, message=FALSE, error=FALSE)
library(patchwork)
library(tidyverse)
library(bayestestR)
library(kableExtra)
library(knitr)
#load("regional-fits.Rdata") # created using save(reg.fit, file="regional-fits.Rdata")
load("report-africa-subregional-fit_31-Mar-2022.Rdata")
source("fitting-functions2.R")  # extraction and other functions

B. Last 5 years regional trend

# Estimate the posterior probability that the trend in the last five years is negative.
# Extract the posterior sample for the marginal mean (unweighted) for the last five years:
# yearP.eff.post.long <- extract.posterior(all.fit, "^YearP.est.MM\\[", index.type="year" )

for(i in c(1,2,3,4)) {

yearP.eff.post.long <- extract.posterior(reg.fit[[i]], "^YearP.est.MM\\[", index.type="year" )
region.name <- reg.fit[[i]]$subregion.index[1,1]

# select the last five years
last.five.years <- sort(unique(yearP.eff.post.long$year),decreasing=TRUE)[1:5]

yearP.eff.post.long <- yearP.eff.post.long[ yearP.eff.post.long$year %in% last.five.years, ]

# compute the slope in pike in the last five years for each sample from the posterior
yearP.eff.slope <- plyr::ddply(yearP.eff.post.long, "sim", function(x){
    fit <- lm( value ~ year, data=x)
    slope <- coef(fit)[2]
    slope.neg <- slope < 0
    data.frame(intercept=coef(fit)[1], slope=slope, slope.neg=slope.neg)
})

# plot estimates and lines 
gg <- ggplot(yearP.eff.post.long, aes(x= year, y=value)) + geom_point(position=position_jitter(w=0.2), alpha=.2) 
gg <- gg + geom_abline(data=yearP.eff.slope, aes(intercept=intercept, slope=slope, x=NULL, y=NULL), alpha=0.01)
gg <- gg + ggtitle(region.name)

# Plot distribution of the post. belief slope with HDI bounds

HDI95 <- ci(yearP.eff.slope$slope,method = "HDI", ci = .95)
HDI95.ci <- data.frame(x=unlist(HDI95)[2:3], y=0)

post.belief.slope.negative <- formatC(mean(yearP.eff.slope$slope.neg), digits=2, format='f')

slope.post <- ggplot(data=yearP.eff.slope, aes(x=slope, y=..density..)) +
              geom_histogram(alpha=0.2) + geom_vline(xintercept=0, color="red") +  
              ggtitle("Posterior distribution of slope in \n fitted yearly PIKE in \n last five years",
              subtitle="Based on unweighted marginal mean\n-MM.p.uw") +
              annotate("text", label=paste("Post belief that slope is < 0 is ", post.belief.slope.negative, sep=""),
              x=Inf, y=Inf, hjust=1, vjust=1) +
              xlab("Slope in estimated PIKE in last five years") + ylab("Density")

slope.post <- slope.post + geom_vline(xintercept = HDI95.ci[,1], lty =2) 
slope.post <- slope.post + geom_text(data=HDI95.ci, aes(x=x, y=y, label=round(x,3))) 
 
# put the two plot together
print(gg + slope.post)



# Print results in table form 

HDI95.tbl.ci <-  yearP.eff.slope %>% 
                 summarize(region = region.name,
                   years = paste(range(last.five.years), collapse = "-"),
                   mean.slope = mean(slope),
                   HDI.lo =  ci(slope ,method = "HDI", ci = 0.95 )[[2]],
                   HDI.hi =  ci(slope ,method = "HDI", ci = 0.95 )[[3]],
                   HDI.ci  = 0.95,
                   HDI.flag = ifelse((HDI.lo <= 0) & (HDI.hi >= 0),
                                        "HDI includes zero", "HDI excludes zero"),
                   post.belief.slope.negative = mean(slope.neg))  %>%
                   ungroup()


print(kable(HDI95.tbl.ci, digits = 3) %>%  
      kable_styling(bootstrap_options = c("striped", "hover")))

cat("\n")
}
region years mean.slope HDI.lo HDI.hi HDI.ci HDI.flag post.belief.slope.negative
Southern Africa 2017-2021 -0.061 -0.086 -0.036 0.95 HDI excludes zero 1
region years mean.slope HDI.lo HDI.hi HDI.ci HDI.flag post.belief.slope.negative
West Africa 2017-2021 -0.02 -0.087 0.044 0.95 HDI includes zero 0.709
region years mean.slope HDI.lo HDI.hi HDI.ci HDI.flag post.belief.slope.negative
Eastern Africa 2017-2021 -0.032 -0.056 -0.009 0.95 HDI excludes zero 0.996
region years mean.slope HDI.lo HDI.hi HDI.ci HDI.flag post.belief.slope.negative
Central Africa 2017-2021 -0.055 -0.093 -0.015 0.95 HDI excludes zero 0.997